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I present a novel algorithm for reconstructing the Wigner 
function from homodyne statistics. The proposed method, 
based on maximum-likelihood estimation, is capable of com- 
pensating for detection losses in a numerically stable way. 

PACS Number(s): 42.50.Dv, 03.65.Bz 



An intriguing aspect of quantum mechanics is the intri- 
cate description of the state of a physical system. There- 
fore, a great deal of interest has been attracted by the re- 
cent experimental demonstration that the quantum state 
of a simple system, namely a single light mode, can be 
completely characterized in a feasible scheme [[y. The 
measurement was based on an observation that marginals 
of the Wigner function, which contains complete informa- 
tion on the quantum state, can be collected with the help 
of a balanced homodyne detector Q . The Wigner func- 
tion was reconstructed from the statistics of homodyne 
events using a standard filtered back-projection algo- 
rithm developed in image processing. This seminal work 
initiated an extensive research in the field of quantum 
state measurement, which has brought beautiful demon- 
strations of quantum phenomena as well as deeper un- 
derstanding of the foundations of quantum theory || . 

The Wigner function is particularly well suited for vi- 
sualization of quantum states, as it represents pictori- 
ally quantum coherence in the form of nonclassical phase 
space structures. The purpose of this Communication is 
to present a novel method for reconstructing the Wigner 
function from homodyne statistics, tts essential advan- 
tage compared to the standard back-projection algorithm 
is the capability of compensating, in a numerically stable 
way, for non-unit efficiency of the homodyne detector, 
tmperfect detection is well known to have a deleterious 
effect on nonclassical features of the Wigner function, 
such as negativities and oscillatory interference patterns 
Furthermore, an attempt to incorporate compensa- 
tion into the standard linear reconstruction scheme fails 
due to rapidly exploding statistical errors and numer- 
ical instabilities [^-^). In the present Communication I 
show that these difficulties can be effectively overcome by 
taking into account a priori constraints imposed by the 
quantum mechanical form of the Wigner function. This 
new method is derived using the maximum-likelihood es- 
timation ||, and constraints are formulated in a way 
which provides a convenient algorithm for reconstructing 
the Wigner function from realistic, finite and imperfect 
homodyne data. 

Let us start the considerations from tracing the stan- 



dard route from raw statistics of homodyne events to the 
Wigner function. The quantum mechanical probability 
distribution of measuring the quadrature x in a single run 
of the homodyne setup with the local oscillator phase 9 
is given by the expectation value of the following positive 
operator- valued measure EQ|: 
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where rj is the efficiency of the homodyne detector, and 
xg is the quadrature operator. In the limit of perfect 
detection, the above expression reduces to 



H(x;6) 



8(x - xg), 
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i.e., the measured statistics become distributions of 
quadrature operators, which are one-dimensional projec- 
tions of the Wigner function. This projection relation 
can be analytically inverted, which yields an expression 
for the Wigner function as the inverse Radon transform 
pT~| of the family of quadrature distributions. 

In standard optical homodyne tomography, regular- 
ized inverse Radon transform is applied to frequency 
histograms of homodyne events. In other words, fre- 
quency histograms are regarded as estimates for quan- 
tum mechanical quadrature distributions, and statisti- 
cal uncertainty in the reconstruction process is governed 
by a simple propagation law. This treatment of statis- 
tical noise has crucial consequences, when efficiency of 
the homodyne detector is less than one. In this case, 
application of the inverse Radon transform to a family 
of distributions given by Eq. (|l|) yields a generalized, s- 
ordered quasidistribution function with the ordering pa- 
rameter s = — (1 — 77) /?7 Q. This object is related to 
the Wigner function by a convolution with a Gaussian 
function. However, the inverse relation is practically use- 
less in numerical calculations: linear deconvolution is an 
ill-posed problem, which would enormously amplify sta- 
tistical noise inevitably present in experimental data [|| . 
The source of this difficulty lies in the estimation proce- 
dure applied in standard optical homodyne tomography. 
This procedure is performed at the level of experimen- 
tal frequency histograms, which are inserted in place of 
quantum mechanical distributions. In such an approach, 
statistical noise is regarded as necessary evil, whose im- 
pact can be quantified, but cannot be reduced. 

Fortunately, this conclusion becomes invalid if we 
adopt a more careful treatment of experimental data. 
In fact, there is a powerful tool for improving the per- 
formance of the reconstruction procedure: it is a priori 
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knowledge about constraints satisfied by quantities to be 
estimated. Use of appropriate estimation methodology 
allows one to incorporate this information in the recon- 
struction scheme. Such an approach can lead to sub- 
stantial reduction of the statistical uncertainty, as in this 
case estimates are picked only from a priori restricted 
physically sensible region. Specifically, the quantum me- 
chanical definition of the Wigner function results in cer- 
tain constraints on its values. The practical problem is 
how to translate this a priori knowledge into an efficient 
numerical algorithm. In the following, I will show that 
this goal can be effectively achieved for optical homodyne 
tomography. 

We will start from the observation that the Wigner 
function at a phase space point (q,p) is given by the 
expectation value of the operator W(q, p) , which can be 
decomposed into the diagonal form |12| : 

i 00 

W{q,p) = -Y j {-l) n Q n [q,p) (3) 

n=0 

where g n {q,p) are projections on displaced Fock states: 

g n (q,p) = D(q,p)\n)(n\D\q,p). (4) 

Here D(q,p) denotes the displacement operator. Quan- 
tum expectation values of g n (q,p) satisfy obvious condi- 
tions: 

00 

(Qn(q,p))>0, Y,{§ n (q,p)) = l (5) 

n=0 

The decomposition defined in Eq. (J3J) suggests the fol- 
lowing two-step reconstruction scheme: for a given phase 
space point (q,p) find estimates for the family of observ- 
ables g n (q,p), taking into account constraints given by 
Eq. (||). Using these estimates, compute the value of the 
Wigner function according to Eq. (||). What makes this 
scheme more robust to statistical noise compared to stan- 
dard optical homodyne tomography, is that estimates of 
(Qn{Q,p)) are a priori restricted to the physically sensible 
region. In contrast, estimates for positive definite observ- 
ables, such as Fock state projections, obtained via stan- 
dard tomographic technique of pattern functions, may, 
and often do take unphysical negative values J?],|| . These 
artifacts become particularly strong when compensation 
for detection losses is built into the pattern functions. 

The positive definite estimates for ((3 n (<?,£>)) will be 
found using the maximum-likelihood approach. The esti- 
mation procedure is motivated by the one-to-one relation 
linking the photon distribution with the phase-averaged 
homodyne statistics jL3|. This relation can be expressed 
in the operator form as 

— / d0H(y;0) = J2My)\n)(n\, (6) 

where 



describes contribution generated by the occupation of the 
nth Fock state. Here Hk(y) is the fcth Hermite polyno- 
mial. The relation given by Eq. (^J), along with statistical 
characterization of experimental data, has been shown to 
provide an algorithm for reconstructing the photon dis- 
tribution with positivity constraints |Q. This approach 
will be now generalized to the estimation of projections 
on arbitrarily displaced Fock states. 

For this purpose let us apply the coherent displacement 
transformation to both the sides of Eq. @, and denote 
the transformed left hand side by T(y; q,p). Making use 
of the identity 

D{q, p)xgD^ (q, p) = xg — q sin 9 — p sin 9 (8) 

yields 

1 [* 

T(y;q,p) = — / d0 H(y + ^/rjq cos 9 + y/rjpsind; 9) 

OO 

= Y^A n {y)Q n {q,p). (9) 

n=0 

Thus, T(y;q,p) as a function of y describes a proba- 
bility distribution obtained by integrating the family of 
shifted homodyne statistics over the local oscillator phase 
9. Eq. (H) shows that T(y; q,p) is uniquely defined by the 
set of projections on displaced Fock states g n (q,p) and 
vice versa: all the quantum expectation values (g n (q,p)} 
can be recovered from (T(y; q,p)}- For 77 = 1, the trans- 
formation of homodyne statistics in Eq. (^|) has a simple 
pictorial interpretation in phase space: moving the phase 
space origin to the point (q,p) corresponds to displacing 
the distribution of the quadrature xg by qcos9 +psiri9. 

Let us now look at the derived relation from the per- 
spective of raw experimental results. Data collected in N 
runs of the homodyne experiment consist of pairs (a^, 9i) 
specifying the outcome Xi and the phase 9i for an ith 
run of the setup fl5|| . For a given phase space point (q,p) 
these data can be converted to the form 

Vi = Xi - ^/rjq cos 0, - ^/rjp sin 4 . (10) 

Statistics of these outcomes is governed by the expecta- 
tion value of T(y;q,p) over the quantum state that is 
measured. This probability distribution (T(y; q,p)) con- 
tains complete information on the set of (g n {q,p))- How- 
ever, in a real experiment the distribution (T(y; q,p)} is 
not known perfectly. The only information we have in 
hand is a set of outcomes characterized by this distribu- 
tion, and a priori knowledge that {g n {q,p)) generating it 
are positive definite and sum up to one. From this in- 
complete information we need to infer estimates for pro- 
jections on displaced Fock states, which we will denote 
for short by g n . The solution given to this problem by 
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the maximum-likelihood methodology is to pick the set 
of {g n } for which it was the most likely to obtain the 
actual result of the series of measurements. Mathemati- 
cally, this is done by maximization of the log-likelihood 
function |B| 



i \ n / n 



(ii) 



where the outcomes yi are treated as fixed parameters for 
a given phase space point (q,p). In the above formula, the 
method of Lagrange multipliers has been used to include 
the constraint Y] g n = 1. 

Thus we have eventually arrived at the following com- 
putational recipe for reconstructing the Wigner func- 
tion from raw homodyne outcomes: for a given phase 
space point (q,p), convert experimental data according 
to Eq. ( |l0| ) and find the maximum of the corresponding 
log-likelihood function £({g n }; {y^}) over the manifold 
defined by Eq. (^). Finally, evaluate the Wigner func- 
tion as 



(12) 



The nontrivial step in the above scheme is the multidi- 
mensional constrained optimization necessary to maxi- 
mize the log-likelihood function. In solving this problem, 
it is useful to note its specific form: the statistics of {y^} 
depends linearly on {g n } which has all the properties of 
a probability distribution as a function of n. This makes 
our task a special case of linear inverse problems with 
positivity constraints [lq|. An effective tool in solving 
this class of inverse problems is the so-called expectation- 
maximization (EM) algorithm jl7|]. Its principle of op- 
eration can be understood by considering the necessary 
condition for the maximum of the log-likelihood func- 
tion £({g n }; {yt})- For each m, the partial derivative 
dC/dg m must vanish, unless the maximum is located on 
the boundary of the allowed region for which g m = 0. 
These two possibilities can be written jointly as 



Qn 



dC 

dg m 



0. 



This condition can be rearranged to the form 

_ J_ \p A m (yj)g m 
0m ~ N^fY, n MVi)Qn 



(13) 



(14) 



for all ms, which shows that the maximum likelihood 
estimate for {g n } is a fixed point of the nonlinear trans- 
formation defined by the right hand side of Eq. (|l4|) . The 
idea of the EM algorithm is simply to iterate this trans- 
formation. When sufficient mathematical conditions 
are fulfilled, this procedure converges to the maximum- 
likelihood solution Jl6| . 

I illustrate the presented method with the reconstruc- 
tion of the Wigner function from Monte Carlo simulated 
imperfect homodyne statistics for a Schrodinger cat state 



y/2 - 2exp(-2|af) 



(l«) 



a)) 



(15) 



with a — 2i. The computer generated data consisted 
of 10 5 events simulated for each of 64 phases uniformly 
spaced between and w. The detector efficiency was 
assumed to equal rj = 90%. The effect of imperfect de- 
tection can be observed in the histogram of homodyne 
events for the phase 9 = 0, depicted in Fig. The visi- 
bility of interference fringes in the homodyne statistics is 
substantially smaller than in the corresponding quadra- 
ture distribution. This blurring is a result of detection 
losses, and it analogously affects the quasidistribution 
function reconstructed by means of inverse Radon trans- 
form, decreasing the magnitude of the oscillatory pattern 
in phase space. 

Fig. H shows the Wigner function reconstructed along 
the position axis q via 10 4 iterations of the EM algorithm 
at each point, starting from a flat distribution of {g n } for 
< n < 39. For a quantitative comparison, the recon- 
structed values are plotted together with the true Wigner 
function evaluated for the state \^f). It is clearly seen that 
virtually full magnitude of the oscillatory pattern is re- 
covered despite detection losses. In contrast, the dashed 
line in Fig. || depicts the quasidistribution function char- 
acterized by the ordering parameter s = — (1 — 77) / 77. This 
function would have been obtained using the standard 
linear back-projection algorithm from homodyne statis- 
tics in the limit of infinite number of measurements. 

In conclusion, the reconstruction algorithm presented 
in this Communication demonstrates that application 
of appropriate estimation methodology can substantially 
improve performance of quantum homodyne tomography. 
The maximum-likelihood approach applied in this paper 
provides an algorithm that is entirely free from singu- 
larities appearing in the standard linear reconstruction 
scheme. Of course, incorporation of a priori constraints 
does not automatically cancel all the effects of imper- 
fect detection. For a fixed number of measurements, the 
quality of the maximum- likelihood estimate worsens with 
decreasing detector efficiency 77. What is important, how- 
ever, is that in any case the obtained estimate for {g n } 
remains between physical bounds, which allows one to 
evaluate safely the sum in Eq. ([12]) . Quantitative discus- 
sion of statistical uncertainty in the proposed algorithm 
will be a subject of a separate publication. 

It is noteworthy that the EM algorithm is a well known 
tool in classical tomography and image restoration p8[ . 
However, these applications make essential use of the fact 
that the multidimensional object to be reconstructed is 
positive definite. This is not the case of the Wigner func- 
tion, whose negativities are a signature of nonclassical 
properties. Therefore, classical reconstruction algorithms 
with a priori constraints usually cannot be directly used 
in quantum state measurement, and there is an appar- 
ent need to develop novel tools for quantum tomography, 
such as that presented in this paper. 
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FIG. 1. Histogram of Monte Carlo simulated homodyne 
events for the local oscillator phase 8 = 0. The modulation 
depth of the interference pattern is decreased by imperfect 
detection. 
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FIG. 2. A comparison of the reconstructed Wigner func- 
tion (•) with its analytical form (solid line). The dashed line 
represents the quasidistribution function rj~ 1 W{r]~ 1 ^ 2 q,0\ 
— (1 — *?)/??) that is related via inverse Radon transform to 
blurred quadrature distributions corresponding to efficiency 
V = 90% §. 



